﻿namespace EnergyGeneratorSim
{
    using System;
    using System.Collections.Generic;
    using System.IO;

    public class Magnetkonverter1
    {
        private const int AnzPmax = 10000;  // Anzahl der Zeitschritte zur Lsg. der Dgl.

        private long n;  //Windungszahl der Spule
        private double A; // .Querschnittsfläche der Spule
        private double Bo; //Magnetfeld (Amplitude) des Dauermagneten
        private double ls; // Länge der zylindrischen Spule
        private double di; // Durchmesser des Spulenkörpers
        private double Dd; // Drahtdurchmesser
        private double rm; // Radius des Dauermagneten
        private double L; // Induktivität der Spule
        private double C; // Kapazität des Kondensators
        private double R; // Ohm`scher Widerstand des Spulendrahtes
        private double rho; // Spez. Widerstand von Kupfer, je nach Temperatur, Kohlrausch,T193

        private double[] phi = new double[AnzPmax];   // Drehwinkel und dessen Ableitungen
        private double[] phip = new double[AnzPmax];
        private double[] phipp = new double[AnzPmax];

        private double[] Q = new double[AnzPmax]; // Ladung und deren Ableitungen
        private double[] Qp = new double[AnzPmax];
        private double[] Qpp = new double[AnzPmax];

        private long AnzP; // Anzahl der tatsächlich berechneten Zeit-Schritte
        private double dt;  // Dauer der Zeitschritte zur Lsg. der Dgl.
        private int Abstd; // Jeder wievielte Punkte soll geplottet werden
        private double omo; // Kreis-Eigenfrequenz des elektrischen Schwingkreises
        private double T; // Schwingungsdauer des elektrischen Schwingkreises
        private double UC; // Kondensatorspannung und Spulenspannung
        private double rhom; // Dichte des Magnetmaterials
        private double dm; // Dicke des zylindrischen Dauermagneten
        private double mt; // Träge Masse des zylindrischen Dauermagneten
        private double J; // Trägheitsmoment des zylindrischen Dauermagneten

        private double[] K0 = new double[AnzPmax];    // Kontroll-Felder für Plot-Zwecke
        private double[] K1 = new double[AnzPmax];
        private double[] K2 = new double[AnzPmax];
        private double[] K3 = new double[AnzPmax];
        private double[] K4 = new double[AnzPmax];
        private double[] K5 = new double[AnzPmax];

        private double EmA, EmE, siA, siE; // Energie: Mittelwerte und Sigma "Anfang" und "Ende"
        private double delE, sigdelE; // Veränderung der Energie-Mittelwerte "Anfang" zu "Ende"
        private double UmAn; // Startwert: Umdrehungen pro Sekunde bei Anlaufen des Rotors
        private double Eent; // Entnommene Energie, elektrisch
        private double Rlast; // Ohm'scher Lastwiderstand

        private void Wait()
        {
            Console.Write("<W>");
            var input = Console.ReadKey();
            Console.Write(input.KeyChar);

            if (input.KeyChar == 'e')
            {
                Environment.Exit(0);
            }
        }

        private void ExcelAusgabe(string fileName, int spalten, IList<double> KA, IList<double> KB, IList<double> KC, IList<double> KD, IList<double> KE, IList<double> KF, IList<double> KG, IList<double> KH, IList<double> KI, IList<double> KJ, IList<double> KK, IList<double> KL)
        {
            using (var stream = File.OpenWrite(fileName))
            using (var fout = new StreamWriter(stream))
            {
                for (int lv = 0; lv < AnzP; lv++)
                {
                    if (lv % Abstd != 0)
                    {
                        continue;
                    }

                    for (int j = 0; j < spalten; j++)
                    {
                        double wert = double.NaN;

                        if (j == 0) wert = KA[lv];
                        if (j == 1) wert = KB[lv];
                        if (j == 3) wert = KC[lv];
                        if (j == 4) wert = KD[lv];
                        if (j == 5) wert = KE[lv];
                        if (j == 6) wert = KF[lv];
                        if (j == 7) wert = KG[lv];
                        if (j == 8) wert = KH[lv];
                        if (j == 9) wert = KI[lv];
                        if (j == 10) wert = KJ[lv];
                        if (j == 11) wert = KK[lv];
                        if (j == 12) wert = KL[lv];

                        fout.Write(wert.ToString().Replace('.', ','));
                        fout.Write('\t');   // Daten-Trennung, Tabulator
                    }

                    fout.WriteLine();   // Zeilen-Trennung
                }
            }
        }

        public void Calculate()
        {
            // Initialisierung - Vorgabe der Werte: Wir arbeiten in SI-Einheiten
            Console.WriteLine("Raumenergie-Konverter mit Rotation.");
            Console.WriteLine("Lichtgeschwindigkeit c = " + Physics.Lichtgesch + " m/s");

            // Spule, Magnet, Kondensator:
            n = 1600; // Windungszahl der Spule
            di = 0.09; // Spulenkörper-Durchmesser
            Dd = 0.0010; // Drahtdurchmesser
            Bo = 0.700; /*Tesla*/ // Magnetfeld (Amplitude) des Dauermagneten
            ls = 0.01; /*Meter*/ // Länge des zylindrischen Spulenkörpers
            C = 0.23E-6; /*Farad*/ // Kapazität des Kondensators
            rm = 0.039; /*Meter*/ // Radius des zylindrischen Dauermagneten
            dm = 0.01; /*Meter*/ // Dicke des zylindrischen Dauermagneten
            rhom = 7.8E3; // Dichte des Magnet-Materials, Eisen, Kohlrausch Bd.3

            // Abgeleitete Parameter, keine Eingabe möglich:
            A = di * di; /*Meter * Meter*/ // Querschnittsfläche der Spule
            L = Physics.Muo * A * n * n / ls; // Induktivität der Spule
            omo = 1 / Math.Sqrt(L * C); // Kreis-Eigenfrequenz des elektrischen Schwingkreises
            T = 2 * Math.PI / omo; //Schwingungsdauer des elektrischen Schwingkreises
            rho = 1.7E-8; /*Ohm*m*/ // Spez. Widerstand von Kupfer, je nach Temperatur, Kohlrausch,T193
            R = rho * (2 * Math.PI * di * n) / (Math.PI * (Dd / 2) * (Dd / 2)); /*Ohm*/ // Ohm`scher Widerstand des Spulendrahtes

            // Sonstige:
            UmAn = 100; // Startwert: Umdrehungen pro Sekunde bei Anlaufen des Rotors
            Rlast = 28; // Ohm'scher Lastwiderstand
            AnzP = AnzPmax; // Anzahl der Zeitschritte insgesamt
            dt = 0.0001; /*sec.*/ // Größe der Zeitschritte
            Abstd = 1; // Jeder wievielte Punkte soll geplottet werden
            mt = Math.PI * rm * rm * dm * rhom; // Träge Masse des zylindrischen Dauermagneten
            J = 1.0 / 2.0 * mt * rm * rm; // Trägheitsmoment des zylindrischen Dauermagneten

            // Anzeige der Werte:
            Console.WriteLine("Induktivitaet der Luft-Spule: L = " + L + " Henry");
            Console.WriteLine("Eigen-Kreisfreq harmon.el.Osz.: omo = " + omo + " Hz => T = " + T + "sec.");
            Console.WriteLine("Laenge des Spulendrahts: " + (2 * Math.PI * di * n) + " m");
            Console.WriteLine("Ohm`scher Widerstand des Spulendrahts: R = " + R + " Ohm");
            Console.WriteLine("Traege Masse des zylindrischen Dauermagneten: mt = " + mt + " kg");
            Console.WriteLine("Traegheitsmoment des Dauermagneten: J = " + J + " kg*m^2");
            Console.WriteLine("Gesamtdauer der Betrachtung: " + AnzP * dt + " sec.");

            // Hier beginnt das Rechenprogramm.
            Console.WriteLine("Mechanisch und elektrisch gekoppelte Schwingung.");
            UC = 0; /*Volt*/
            Q[0] = C * UC;
            Qpp[0] = 0;
            Qp[0] = 0; // Elektrische Startwerte
            phi[0] = 0;
            phip[0] = UmAn * 2 * Math.PI;
            phipp[0] = 0; // Mechanische Startwerte
            Eent = 0; // Reset für: Entnommene Energie, elektrisch
            K0[0] = 0;
            K1[0] = 1.0 / 2.0 * L * Math.Sqrt(Qp[0]); // Spulen-Energie
            K2[0] = 1.0 / 2.0 * C * Physics.Sqr(Q[0] / C); // Kondensator-Energie
            K3[0] = 1.0 / 2.0 * J * Physics.Sqr(phip[0]); // Rotations-Energie
            K4[0] = K1[0] + K2[0] + K3[0]; // Gesamt-Energie
            K5[0] = 0;

            for (int i = 1; i < AnzP; i++)
            {
                Qpp[i] = -1.0 / L / C * Q[i - 1] - (R + Rlast) / 2 / L * Qp[i - 1]; // nach *5 von S.6
                Qpp[i] = Qpp[i] + n * Bo * A * Math.Sin(phi[i - 1]) * phip[i - 1] / L; // Induzierte Spannung in Spule bringen.
                Qp[i] = Qp[i - 1] + (Qpp[i] - R / 2 / L * Qp[i - 1]) * dt; // nach *3 & *4 von S.6
                Q[i] = Q[i - 1] + Qp[i] * dt;
                phipp[i] = -Bo * n * Qp[i] * A / J * Math.Sin(phi[i - 1]); // Mechanisches Drehmoment, x-Komponente
                phip[i] = phip[i - 1] + phipp[i] * dt;
                phi[i] = phi[i - 1] + phip[i] * dt;
                K0[i] = 0;
                K1[i] = 1.0 / 2.0 * L * Physics.Sqr(Qp[i]); // Spulen-Energie
                K2[i] = 1.0 / 2.0 * C * Physics.Sqr(Q[i] / C); // Kondensator-Energie
                K3[i] = 1.0 / 2.0 * J * Physics.Sqr(phip[i]); // Rotations-Energie
                K4[i] = K1[i] + K2[i] + K3[i]; // {Gesamt-Energie
                K5[i] = Rlast * Physics.Sqr(Qp[i]); // Am Lastwiderstand entnommene Leistung
                Eent = Eent + K5[i] * dt; // Am Lastwiderstand entnommene Energie
            }

            // Gesamt-Energie-Bilanz und Anzeige:
            EmA = 0;
            EmE = 0;
            siA = 0;
            siE = 0;

            for (int i = 0; i < 10; i++)
                EmA = EmA + K4[i] / 10; // Mittelwert zu Beginn
            for (int i = (int)AnzP - 9; i < AnzP; i++)
                EmE = EmE + K4[i] / 10; // Mittelwert am Ende
            for (int i = 0; i < 10; i++)
                siA = siA + Physics.Sqr(EmA - K4[i]); // Varianz zu Beginn
            for (int i = (int)AnzP - 9; i < AnzP; i++)
                siE = siE + Physics.Sqr(EmE - K4[i]); // Varianz am Ende

            siA = Math.Sqrt(siA) / 10;
            siE = Math.Sqrt(siE) / 10; // Standardabweichungen

            Console.WriteLine("Energie-Werte: E_Anfang = (" + EmA + " +/- " + siA + ") Joule");
            Console.WriteLine(" E_Ende = (" + EmE + " +/- " + siE + ") Joule");
            delE = EmE - EmA;
            sigdelE = Math.Sqrt(Physics.Sqr(siE) + Physics.Sqr(siA));
            Console.WriteLine("=> Veraenderung: delta_E = (" + delE + " +/- " + sigdelE + ") Joule");
            Console.WriteLine("=> Konvertierte Leistung = (" + delE / (AnzP * dt) + " +/- " + sigdelE / (AnzP * dt) + ") Watt");
            Console.WriteLine("An Rlast entnom.Leistung = " + Eent / (AnzP * dt) + " Watt");
            ExcelAusgabe("test_04.dat", 12, Q, Qp, Qpp, phi, phip, phipp, K0, K1, K2, K3, K4, K5);
            Wait();
            Wait();
        }
    }
}
